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Abstract 

The computation of two Bayesian predictive distributions which are 
discrete mixtures of incomplete beta functions is considered. The num- 
ber of iterations can easily become large for these distributions and thus, 
the accuracy of the result can be questionable. Therefore, existing algo- 
rithms for that class of mixtures are improved by introducing round-off 
error calculation into the stopping rule. A further simple modification is 
proposed to deal with possible underflows that may prevent recurrence to 
work properly. 

Keywords: Predictive distribution; Bayesian approach; Round-off error; Incomplete 
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1 Introduction 



The K -square and K-prime distributions have been introduced in iLecoutre 



. They can be characterized as mixtures of the c lassical noncentral F and 



noncentral t distributions respectively ( Lecoutrel 199°J ). These two distributions 



are involved in the Bayesian predictive approach for planning and monitoring 



experiments (|Lecoutrel 120011 ). In particular, they are useful tools for sample 
size determination, using the predictive distributions of the test statistics and 
of the limits of confidence intervals under standard normal models, assuming a 
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conjugate prior. It must also be noted that they include as particular cases the 
distributions of the square of the sample multiple correlation coefficient and of 
the sample correlation coefficient. The aim of this article is to provide efficient 
algorithms for the calculation of their cumulative distribution functions (cdfs). 
These cdfs can be expressed in terms of infinite series of multiples of incomplete 
beta function ratios, thus adequate for recursive calculations. More precisely, 
both imply the general form 

oo 

(i) 

3=0 

with 

oo 

s = ±l, 0< 5j <1Vj, = 1 

3=0 

and where Hj(x) involves only the incomplete beta function. 

Dealing with a related problem, the Applied Statistics algorithm AS 278 



developed for the psi-square distribution (jLecoutre. Guigues and Poitevineau . 



1992f ) could be adapted to match the present cdfs. However, AS 278 is a 



Metho d 1 recursive algorithm, in the terms of iBenton and Krishnamoorthv 



( 20031) : accumulation is simply done from index until a convergence crite- 



rion is met. In some cases (especially when the noncentrality parameter of the 
distribution is large), it can lead to an exceedingly large number of itera tions 



and c onsequently to unacceptable execution time and loss of precision. iFrick 
(1990) proposed an improvement that consists in starting iterations at an index 
such that the resulting truncation error is negligible, but this does not solve the 
problem. 



Yet, the present cdfs are of the general class considered bv lBenton and Krishnamoorthv 



(2003) and, as such, are good candidates for what they called Method 2 class of 



algorithms. Essentially, this Method 2 is a both backward and forward recur- 
sive algorithm where the starting index for iterations, say fc, is chosen so that 
gk is a maximum, which reduces the above mentioned problems. Nevertheless, 
although smaller than with Method 1, the number of iterations can still remain 
important as soon as parameters increase. Thus, when a relatively high degree 
of accuracy is required, the problem of round-off errors cannot be neglected. 

Therefore, we present in the next two sections a Method 2 class of algorithms 
that includes round-off error calculations. It is applied here respectively to the 
K-square and K-prime cdfs, but is of general use as far as the general form (TT]) is 
concerned. CPU times are presented in section[4] along with a few illustrations, 
and some examples of applications of these cdfs are given in section [5] In 
section [5] we discuss some remaining problems and propose, in some cases, a 
simple modification which leads to an algorithm that is intermediate between 
Method 1 and Method 2. Section [7] is devoted to some concluding remarks. 
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2 K-square distribution 



Techn ical characterizations of the K-square distribution can be found in lLecoutre 
(1999). This distribution is written K 2 r (a 2 ) where p, q, r are degrees of free- 
dom parameters and a 2 is a noncentrality parameter. 

Particular cases of the K-square distribution are: 
a = : K 2 r (0) = F PtT (usual F distribution), 
q = oo : Kp oar (a 2 ) = Fp r (a 2 ) (noncentral F distribution), 
r = oo : K 2 q oo (a 2 ) = A 2 q (a 2 ) (lambda-square or alternate chi-square distri- 
bution), 

q = oo,r = oo : K 2 tOO OC (a 2 ) = (l/p)xp{a 2 ) (noncentral chi-square distribu- 
tion). 

For the cdf, s = 1 in ([1]) and we simply have 

oo 

Pr(K 2 ^ r (a 2 )<x)=Y,9jH J (x), 
3=0 

with 



93 r(j + l)r(f) \q + a 2 ) \q + a 2 ) (2) 
and 

Hj(x) = I px /( r +px) (| +3, , 0, (3) 
where I z is the incomplete beta function 

Iz{aM = ^+A /V^i -t) 6 - 1 *. 



r(o)r(6) J 

The coefficients gj are the probabilities of obtaining the value j for a vari- 
ate following a negative binomial distribution with parameters q/(q + a 2 ) and 
q/2. The mode is \a 2 (q — 2)/(2 q)\, where [.] denotes the integer part (e.g., see 
Johnson. Kotz and Kemp . 19931 p. 209), hence the starting index for iterations. 



From this, it is clear that the number of iterations heavily depends on a 2 . 

The recurrence relations for the cdf are straightforward. For the H^s (the 
incomplete beta function) we have 



T(p/2 + r/2 + j) ( px 

Hj + i - tij - 



p/2+j , \ r/2 



Hi -i — Hi 



T{p/2 + j + l)r(r/2) \r+pxj \r + px 
r(p/2 + r/2 + j - 1) ( px \ 



r(p/2 + j)T(r/2) \r+px) \r + px 

and for the gj coefficients 

q/2 + j _a 

j + 1 q + a 



p/2+j-l / r N r/2 



j 9 + a 2 
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Let A and 6 denote the absolute and the relative error respectively. The 
absolute error for an individual term of the scries is 

Now, noting k the starting index of the computations, the forward and backward 
recurrences for gj are respectively of the form 



9k+j = gk+j-lCk+j-1 = 9k Y[ c k+i an d 9k-j = 9k-j+l/c k -j = 9k JJ , 

1=0 1=1 Ck - 1 

so that 

i-i j 
$9k+ = Sg k + ^ Sc k+i and 6g k -j = 5g k + ^ 5c k -i. 

i=0 i=l 

If we assume that the relative errors on the coefficients Cj are constant, say equal 
to e (e.g., we can assume that all Cj's are calculated with a maximal precision 
of n decimal digits so that e < ilO~™ +1 ), we obtain 

Sg k ±j = Sg k +je hence Ag k±j = (5g k + je)g k±j . 

For the terms Hj (x) , the recurrence involves a sum 

i-i 

H k+ j = H k+j -i - d k+ j-i = H k — d k +i, 

i=0 

j 

Hk-j = H k _j + \ + d k -j = H k + d k -i, 

i=i 

then, 

j-i 3 
AH k+j = AH k + ^Adh+i and AH k -j = AH k + ^Ad fc - 4 - 

*=0 i=l 

The coefficients d k ±i contain gamma functions which can themselves be calcu- 
lated by recurrence, just as for the g^s. Therefore, with the same assumptions 
as for the coefficients gj, we have 

Ad k±j = (8d k +je)d k ±j. 

Consequently, the round-off error (E c ) of a calculation involving N iterations 
(both backward and forward) becomes 

N min(JV,fe) 

E c = A(g k H k ) + J2 A (9k+jH k+j )+ &{9k-jH k -j). (4) 
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Skipping tedious but elementary calculations, it gives 



E c = (Sgk + SH k )g k H k + 

N 

22{(SH k + 6d k )gk+jHk + (Sgk + je - Sd k )g k+j H k+:j + 
3=1 

3-1 

e 9k+j^2id k+l } + 

i=0 
min(JV,fc) 

{(SHk - Sd k )g k -jH k + (5g k + je + Sdk )g f , ///. , + 

3=1 

3 

tgk-j^idk-i}- (5) 

z=l 

Now, for the same reason as for the relative errors on the coefficients Cj, we 
can assume Sg k = 8d k — e. Furthermore, H k involves only one calculation of 
the incomplete beta function for which there exist very performing algorithms 
(e.g., AS 63 bv lMaiumder and Bhattachari so that, again, SH k = e is 

a reasonable assumption. Consequently, it reduces finally to 



E r _ = e 



N N N f 3-1 "J 

2H k ^2 9k+] + ^2 i9k+jH k+:j + < gk+j ^2 idk+t \ 

3=0 3=1 3=1 I 1=0 ) 



3=0 3=1 3= 

mxa(N,k) min(JV,fc) 

2 ^2 9k-jHk-j + ^2 39k-jHk-j- 

3=1 3 = 1 

roin(iV,fe) ( j 

Y \ 9k-3 id k-i 
3 = 1 I i=l 



(6) 



Given the parameters, Hj{x) is a decreasing function of j. Thus, when stop- 
ping the calculations at step j, the truncation error (E t ) is bounded by: 

while j < k 



k—j — l oo 

E t < H (x) J2 9r+H k {x) J2 9i 

i=0 i=k+j + l 

k—j — 1 oo 

< H (x) ^2 9i+H Q (x) ^2 9i 

i=0 i=k+j + l 

k+3 

< H (x) 1- J2 Si 

i=k-j 



(7) 



5 



and when j > k 



Et < H k+j {x) 



1 



k+j 

E 

i=0 



Si 



(8) 



is a slight modification of the rule in step 3 in lBenton and Krishnamoorthv 
(|2003f ) who used 1 instead of Hq(x). The relaxation of the stopping rule com- 
pensates for the increased execution time due to one call to the incomplete beta 
function. 

Stopping rule: Stop when E t + E c becomes lower than a predetermined 
absolute error bound or when E c exceeds that error bound, which means that 
the required accuracy cannot be reached. 

For the distribution of the square of the sample multiple correlation coeffi- 
cient (see end of section \Sj , we compared the algorit hm for the K-square cdf, 
called K2CDF, to lBenton and Krishnamoorthvl (|2003l) Algorithm 7.1 (the mode 
of the negative binomial distribution, instead of the mean, was used as the start- 
ing point to ensure the comparability of the two algorithms) . For the examples 
in their Table 1, all results agreed within the 10~ 12 limit that was chosen as 
the maximum absolute error parameter (both algorithms were run in "double 
precision" , i.e. 64-bit words). 



3 K-prime distribution 



Techn ical characterizations of the K-prime distribution can be found in lLecoutre 
(|l999h . This distribution is written K' (a) where q,r are degrees of freedom 



parameters and a is a noncentrality parameter. 

Particular cases of the K-prime distributions are: 
a = : K' q r (0) = t r (usual t distribution), 
q = oo : K'qq r (a) = t' r (a) (noncentral t distribution), 
r = oo : K' q oo (a) = A' (a) (lambda-prime distribution), 
q = oo, r = oo : K'^ ^(a 2 ) = N(a, 1) (normal distribution). 

This cdf has the following properties: 

Fr(K' g<r (a)<x)=Pr(K! riq (x)>a), 

Pr(X; r (-a) < -x) = ?*{KA X ) > a )> 
Pr(Jf; r (a) < 0) = Pr(A' g (a) < 0) = Pr(i, > a). 
Several cases are to be distinguished for the cdf: 

If a > and x < 

Pr(K' q<r {a) < x) = Vv{K' q r (a) < 0) - Pr(a; < K' q r (a) < 0) 



Pr(t q > a) - ^{-lYgjI^/^^ 
3=0 



j + 1 r\ 
2 ' 2 ' 



G 



where 

i r(g±i) / q y / a* \i 
9 > 2r(i±i)r(§) U + «V U + «V ■ 1 ' 

If a > and x > 

Pr(X; r (a) < s) = Pr[K' qr {a) < 0) + Pr(0 < K' q r {a) < x) 



Pr(t 9 > a) + > ,9jlx^/(r+x 2 ) 

3=0 



3 + 1 r 
2 ' 2 



If a < 0, we reduce to the above cases using 

Pr{K' qr {a) <x) = l- Pr{K' qr (-a) < -x). 
If a = 0, we simply have 

Pr(i^ r (0) <x) = Pr(i r < a?). 

Hence, the cdf of the K-prime involves the calculation of the cdf of the usual 
Student's t distribution and a series of the general form JT]). The case where 
a and x are of a different sign is an unfavorable one, since the series is then 
alternate. Therefore, in the algorithm called KPPJMECDF, the even and odd 
terms of the series are accumulated separately in order to minimize the number 
of subtractions. 

The recurrence relations for the incomplete beta function now write 

r( 2±i±i) 



J+ z 3 TV J+3 



r(l±i)r(|) \r + x V \r + x< 

r(2±pi) 



r(^±i)r(§) V^ + . 

and for the gj coefficients 



g + 3 

j + 2 q + a 2 



9j+2 - — — 9j, 



3 Q + a 2 

The starting point for iterations is taken as the mode of the gj's, i.e. k — 
[a 2 (q—2) I q] . Again, a 2 is an important factor regarding the number of iterations. 
The calculation of errors developed for the K-square series directly applies here, 
and the stopping rule is the same. 
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4 Numerical examples and CPU time 



Some numerical examples, also illustrating the speed of the algorithms, are 
presented in Tables 1 and 2. The probabilities presented have been calculated 
with a required accuracy of 10~ 4 . In order to estimate the loss of speed due 
to the calculation of round-off errors, we also computed the cdf using only the 
truncation error in the stopping rule to serve as reference CPU times. In the 
last column of the tables, the time increase is expressed as a percentage of these 
reference CPU times. The programs were compiled with the GNU g95 Fortran 
compiler (GCC 4.0.3, Apr. 19 2006), using "standard real" data type (i.e., 
32-bit words), and CPU time was computed through the Fortran CPU.TIME 
subroutine. The programs were run on an Intel M750 1.86 GHz PC (each 
calculation was computed 20,000 times in order to provide a substantial CPU 
time). 

On the one hand, and as easily predictable from the algorithm, it appears 
that calculation of round-off errors is time consuming. On the other hand, ex- 
amples of its usefulness can be given. For that purpose, keeping the required 
accuracy to 10 -4 , we consider that the same algorithm run in "double pre- 
cision" (64-bit words) with an accuracy parameter set to 10~ 9 provides the 
"exact" value. The absolute difference between this reference value and the 
value returned by the algorithm without round-off error calculations is termed 
"error" in the following (the "exact" value is reported in square brackets). In 
all these cases the algorithm with round-off error calculations rightly returns an 
error message indicating the required accuracy cannot be met. 

For the K-square cdf: 
x = 90, p = 10, q = 15, r = 20, a 2 = 10 3 : error = 1.7 x 10~ 4 [0.4168], 
x = 15, p = 10, q = 20, r = 10 5 , a 2 = 80 : error = 6.0 x 10~ 4 [0.9577], 
x = 9, p = 10, q = 100, r = 10 5 , a 2 = 80 : error = 1.2 x 10~ 2 [0.5259]. 

For the K-prime cdf: 
x = 100, q = 10, r = 20, a = 80 : error = 9.0 x 10~ 4 [0.8101], 
x = 20, q = 10, r = 10 5 , a = 20 : error = 4.9 x 10~ 3 [0.5574], 
x = 20.5, q = 200, r = 10 6 , a = 21: error = 1.5 x 10" 1 [0.3730]. 

All these examples involve the largeness of at least one parameter, precisely 
because it is in such cases that the precision of the result may be suspected. An 
illustrated example for the K-prime cdf is presented in the next section. 

5 Examples of applications 

As an illustration of the use of the K-prime and K-square distributions, consider 
the sample size determination under usual normal models. For instance, a simple 
two-sample experiment is designed to compare a new drug with a placebo. The 
goals of the experiment specify that the new drug is considered as effective if the 
raw difference 5 = /j,d — MP i s more than +3. For this purpose, the investigators 
plan to use a two-sample shifted t test with equal numbers of subjects n in each 
group, in order to test H : S = +3 against the alternative Hi : S > +3. Hence, 
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Tabic 1: Time comparison between Algorithm K2CDF and the same algorithm 
without round-off error calculation for computing Vv{K 2 qr (a 2 ) < x) 20,000 
times (time in second) 



X 


P 


q 


r 


a a 


Pr(^ r (a*) < x) 


CPU time 


time increase 


3 


5 


5 


5 


5 


0.6664 


0.20 


08% 


1 


5 


5 


9 


10 


0.1195 


0.11 


17% 


10 


5 


5 


9 


10 


0.9440 


0.14 


29% 


10 


5 


5 


9 


100 


0.2142 


0.25 


14% 


100 


9 


5 


5 


100 


0.9819 


0.53 


31% 


80 


10 


20 


25 


1000 


0.3015 


1.31 


27% 



Tabic 2: Time comparison between Algorithm KPRIMECDF and the same 
algorithm without round-off error calculation for computing Pr(K' q r (a) < x) 
20,000 times (time in second) 



x q 


r 


a 


Pr(^» < x) 


CPU time 


time increase 


-5 5 


5 


0.5 


0.0007 


0.50 


07% 


5 5 


5 


5 


0.5000 


0.34 


10% 


9 5 


5 


5 


0.8763 


0.55 


25% 


5 5 


5 


10 


0.0872 


0.47 


15% 


9 5 


5 


10 


0.4137 


0.77 


26% 


9 5 


10000 


5 


0.9856 


0.45 


16% 


-15 5 


10 


-50 


0.9918 


4.47 


30% 
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the efficacy of the drug will be assessed if 



d-3 

t — 7== > t q fi.05, 



where d is the observed difference, s is the pooled estimate of the common stan- 
dard deviation a and t q 05 is the 5% upper point of the Student's distribution 
with q = 2n — 2 degrees of freedom. 

Suppose that a conjugate prior distribution has been chosen, such as 5\a ~ 
N(do, (2/no)cr 2 ) and a 2 ~ So(Xg ) _1, For instance, this prior can be the pos- 
terior distribution from a pilot study (starting with a noninformative prior). 
Then, for any given sample size n, the probability of achieving the study target 
can be computed from a K-prime distribution, using the predictive distribution 
of the t test statistic: 

* ~ y/l + n/n K' I *° ) , where t = ^ -=^= . 

Suppose that d — +4.35, so = 2.07, n = 10, hence q n = 18 and to = 
+ 1.458. For instance we find for n = 50 the predictive probability: 

Pr(t > +1.6606) = Pr \k' 1S 98 (lA58y/5/6) > 1.6606/V6 ] = 0.7327. 



In order to get predictive probabilities equal to 0.80 and to 0.90, n = 97 and 
n = 1930 subjects in each group are respectively needed. 

Equivalently, the investigators could compute a 90% confidence interval for 
8 and assess the efficacy of the drug if its lower limit is larger than +3. The 
predictive distribution for this lower limit £ = d — t g ,o.05 Syj2/n also involves a 
K-prime distribution: 

l~ do - W2M> + 2/n K>. qo ( ~ 50 ^ 05 ) . 

■ \y/l + n /nj 

Of course, for any fixed n, we find again the same predictive probabilities. 
This is due to the following fundamental property of the cdf (Lecoutre, 1999): 

Pr [K' q0tq (a)<x)=Pv (k'^x) > a) . 

The K-prime distribution can also be used to make predictive statements 
about the standardized difference d/s in a future sample. In the same situation 
as above (two groups with a same sample size) we have: 



d /2(n + n) , f d / n n 



s V no n qo,q \ so y 2(no + n) 

When q — > oo, this distribution tends to the distribution of the parameter 
(5/cr. Thus, with a very large value of n, it could be used to get a statement 
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about the population standardized difference (as an alternative to the A-prime 
cdf). 

For instance, suppose that oIq/sq — 3 and no = 100. Then, taking n = 
500000, 

Pr{d/s > 2.731804) = 1 - Pr [#198,999998 (21.21108) < 19.31484] = 0.9000. 

But, actually, KPRIMECDF cannot provide a sufficiently accurate answer, 
even when the maximum absolute error parameter is set to 10 -2 , and issues an 
error message, while the algorithm without round-off error calculation returns 
a value (0.92) which is in error by 2 times the required accuracy. 

Concerning the K-square distribution, it can be used for the sample size 
determination in ANOVA designs. For instance, a simple g-sample experiment 
is designed to test the equality of g means. A pilot study has already been 
conducted with g groups of equal sample size no, and a F ratio Fq has been 
obtained (under the usual normal model). Assuming an initial non informa- 
tive prior, the posterior predictive distribution for the F ratio in the planned 
experiment with n subjects in each group is a K-square distribution: 

1 + n/n ^ 2 (9-^ 



F ~ 3-1 K B-l,gno-g,gn-g { 1 + no / n F ° 

Suppose that g = 3, uq = 10, Fq = 3.6 and n = 30. Then, given the first 
results, F is distributed as 2iTf 27 87 (5.4) and the probability of obtaining a 
significant F test at 0.05 level is Pr(F > 3.1013) = 0.7792. In order to get 
predictive probabilities equal to 0.80 and to 0.90, n = 33 and n = 54 subjects 
in each group are respectively needed. 

Other uses of the K-prime and K-square distributions are the computation 
of the cdf of the sampling distributions of correlation coefficients. The cdf of 
the sample coefficient r, involving a sample of n independent observations from 
a bivariate normal population with population coefficient p, is a particular case 
of the K-primc distribution: 



Pr(r < x) = Pr 



The cdf of the square of the sample coefficient R 2 , involving a sample of n inde- 
pendent observations from a p-variate normal population with square multiple 
correlation coefficient p 2 , is a particular case of the K-square distribution: 



Pr(i? 2 < x) = Pr 



A p-l,n-l,n-p \ \ n i Ji I2 I 



1 — p 2 J p — 1 1 



6 Limitations and possible improvements 

Draw backs of Method 1 algorithms (in the terms of iBenton and Krishnamoorthv . 
20031) led to the development of Method 2 algorithms. In Method 1, the itera- 
tions start at index j = which maximizes Hj (x) , while in Method 2 they start 
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at index j — k which maximizes gj. Nevertheless, the latter is not systematically 
better. For instance, it can happen that the initial recurrence increment for the 
Hj's is too small with respect to the machine limit so that a zero is returned 
and recurrence is impossible: e.g., for the K-square cdf, this increment term is 
lower than 1CT 307 when p = 10, q = 20, r = 30, a 2 = 500 and x = 0.1. More 
generally, whenever Hk(x) tends to zero quickly with respect to k, Method 1 
algorithms perform better than Method 2 algorithms, because only the first 
terms of the series (JTJ contribute significantly to the sum. And when Hk(x) 
is still close to Ho(x), Method 2 is quasi optimum (with the same parameters 
as in the preceding example, this is the case when x = 99 : Ho (99) » 1 and 
#22 5 (99) = 0.994). 

Obviously, the best method would be to start iterations at the index (between 
and k) which maximizes the product gjHj(x) and not only one of the terms. 
However, this is not easy to determine in general. A tempting solution, when 
Hk{x) is considered too small, would be to choose the modified index, say k', 
such that Hk> (x) reaches a predetermined value; unfortunately, such an inversion 
of the beta cdf involves an iterative procedure and so is to be discarded on 
grounds of speed efficiency. As an alternative, we propose to simply lower k by 
multiplying it by the argument of the incomplete beta function (px/(px + r) for 
the K-square and x 2 jix 2 + r) for the K- prime). 

For example, for the distribution #10,80,200 (500), when x takes the values 
35, 30, 25, and 22, the number of iterations is always 202 (for a precision of 
10~ 4 ), while when turning to the modified starting index, it drops respectively 
to 155, 146, 136 and 128. 

7 Concluding remarks 

We presented an algorithm for two Bayesian predictive distributions of impor- 
tance for monitoring experiments. This algorithm includes round-off error cal- 
culation and is applicable to any cumulative distribution function that can be 
expressed as a discrete mixture of continuous distributions such that the recur- 
rence relation for the discrete coefficients is multiplicative and the recurrence 
relation for the continuous distribution is additive. However, this kind of er- 
ror calculation (which is only an approximation, of course) is time consuming, 
and when speed is a crucial factor, it has to be introduced only when deemed 
necessary. It will be the case, for example, when the required accuracy is high 
and/or when the number of iterations is large so that the precision of the result 
may be suspected. In this regard, the material used (computer and compiler) 
is of importance, particularly through the variable noted e, the precision of 
an "elementary" recurrence calculation. For instance, two different comput- 
ers/compilers storing variables into words of the same size could have different 
e if they use registers of different size to perform computations. We also con- 
sidered the case where the starting index of iterations is such that recurrence 
is impossible due to underflows. The proposed solution, which is an approach 
to the problem of finding the optimum starting index, is to lower this index by 
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a quantity which is the argument of the incomplete beta function, a choice we 
made on empirical grounds and that is likely to be improved. 

References 

Benton, D., Krishnamoorthy, K., 2003. Computing discrete mixtures of contin- 
uous distributions: noncentral chisquare, noncentral t and the distribution of 
the square of the sample multiple correlation coefficient. Comput. Stat. Data 
An. 43, 249-267. 

Frick, H., 1990. A remark on Algorithm AS 226: Computing noncentral beta 
probabilities. Appl. Statist. 36, 311-312. 

Johnson, N.L., Kotz, S., Kemp, A.W., 1993. Univariate Discrete Distributions. 
Second edition. Wiley, New York. 

Lecoutre, B., 1984. L'Analyse Bayesienne des Comparaisons. Presses Universi- 
taircs dc Lille, Lille. 

Lecoutre, B., 1999. Two usefull distributions for Bayesian predictive procedures 
under normal models. J. Statist. Plann. Inference 79, 93-105. 

Lecoutre, B., 2001. Bayesian predictive procedure for designing and monitor- 
ing experiments. In Bayesian Methods with Applications to Science, Policy 
and Official Statistics, Luxembourg: Office for Official Publications of the 
European Communities, 301-310. 

Lecoutre, B., Guigucs, J.-L., Poitcvineau, J., 1992. Distribution of quadratic 
forms of multivariate Student variables. Appl. Statist. 41, 617-627. 

Majumder, K.L., Bhattacharjee, CP., 1973. The incomplete beta integral. Appl. 
Statist. 22, 409-411. 



13 



